# R script for replicating results in the main text
# Zhai & Garside 2022
# original device: MacPro 13, R 4.1.2
# recommended working directory: Desktop

setwd("~/Desktop") # default wkd
rm(list = ls())

# pkgs --------------------------------------------------------------------

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, magrittr, fixest, texreg, ggthemes)

# helpers -----------------------------------------------------------------

source("R_00_helpers.R")

# graphix -----------------------------------------------------------------

theme.default <- theme_get()
theme_set(theme_minimal(base_size = 14))

# data --------------------------------------------------------------------

df.main <- read.csv("data_vote_main.csv", row.names = 1) %>% 
  group_by(land_name) %>% 
  mutate(flooded_land = max(flooded, na.rm = TRUE)) %>% 
  ungroup() %>% 
  filter(flooded_land == 1) 

# table 1 -----------------------------------------------------------------

xx.h1 <- "+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(date)"
fm.h1 <- "v_green_pct~factor(flooded)*factor(date)"
fm.h1.x <- paste0(fm.h1, xx.h1)
fm2.h1 <- gsub("flooded", "flooded2", fm.h1)
fm2.h1.x <- paste0(fm2.h1, xx.h1)

tab1 <- map(
  list(fm.h1, fm.h1.x, fm2.h1, fm2.h1.x),
  ~fit_twfe(df = df.main, fm = .x, fe = "gen_id", cl = "gen_id"))

texreg::screenreg(tab1, digits = 3)

# figure 2 ----------------------------------------------------------------

df.main.h1a <- 
  df.main %>% 
  mutate(nsevere = as.integer(flooded==1 & severe==0),
         nsevere2 = as.integer(flooded2==1 & severe2==0))

fm.h1a <- "v_green_pct~factor(severe)*factor(date)+factor(nsevere)*factor(date)"
fm.h1a.x <- paste0(fm.h1a, xx.h1)
fm2.h1a <- gsub("severe", "severe2", fm.h1a)
fm2.h1a.x <- paste0(fm2.h1a, xx.h1)

fig2.fits <- map(
  list(fm.h1a, fm.h1a.x, fm2.h1a, fm2.h1a.x),
  ~fit_twfe(df = df.main.h1a, fm = .x, fe = "gen_id", cl = "gen_id"))

fig2.fits.tidy <- 
  fig2.fits %>% 
  map(~{coeftable(.x) %>% as.data.frame}) %>% 
  map(~{rownames_to_column(.x, var = "term") %>% select(term, mean = Estimate, se = `Std. Error`)}) %>% 
  map(~{filter(.x, term %in% c(
    "factor(date)2021-09-26:factor(nsevere)1",
    "factor(date)2021-09-26:factor(nsevere2)1",
    "factor(severe)1:factor(date)2021-09-26",
    "factor(severe2)1:factor(date)2021-09-26",
    "factor(date)2021-09-26:factor(nsevere)1",
    "factor(flooded2)1:factor(date)2021-09-26",
    "factor(date)2021-09-26:flooded2_w")) %>%
      mutate(term = case_when(
        grepl("(flooded2|nsevere2)", term) ~ "Secondary", 
        grepl("(flooded|nsevere)", term) ~ "Primary",
        grepl("severe2", term) ~ "Secondary+",
        grepl("severe", term) ~ "Primary+"
      ))}) %>% 
  bind_rows %>% 
  add_column(model = rep(rep(c("Base", "Full"), each=2), 2)) %>% 
  mutate(lwr = mean-1.96*se, upr = mean+1.96*se) %>% 
  mutate(term2 = case_when(
    grepl("Primary", term) ~ "Primary",
    grepl("Secondary", term) ~ "Secondary"
  )) %>% 
  mutate(term3 = case_when(
    grepl("\\+", term) ~ "High",
    grepl("Primary|Secondary", term) ~ "Low"
  ) %>% factor(levels = c("Low","High")))

ggplot(fig2.fits.tidy) +
  geom_pointrange(aes(x="bhat", y=mean, ymin=lwr, ymax=upr, color=model, shape = term3), position = position_dodge(0.7)) +
  facet_grid(rows = vars(term2)) +
  geom_hline(yintercept = 0, color = "darkred", lty = "dotted") +
  coord_flip() +
  scale_color_stata() +
  theme_minimal(base_size = 14) +
  labs(x = NULL, y = NULL, color = "Model", shape = "Severity", title = "Differential Effect") +
  scale_y_continuous(limits = c(-0.01, 0.06), breaks = scales::pretty_breaks(n=5)) +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank())
ggsave("~/Desktop/main_text_figure2.pdf", width = 7, height = 4) 

# cleanup -----------------------------------------------------------------

pacman::p_unload("all")
rm(list = ls())
